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We propose a method to study the time evolution of Bose 
condensed gases perturbed from an initial thermal equilib- 
rium, based on the Wigner representation of the iV-body den- 
sity operator. We show how to generate a collection of ran- 
dom classical fields sampling the initial Wigner distribution in 
the number conserving Bogoliubov approximation. The fields 
are then evolved with the time dependent Gross-Pitaevskii 
equation. We illustrate the method with the damping of a 
collective excitation of a one-dimensional Bose gas. 
PACS: 03.75.Fi, 05.10.Gg, 42.50.-p 

Since the first experimental demonstrations of Bose- 
Einstein condensation in atomic gases p| , the role played 
by finite temperature effects in the physics of Bose con- 
densed alkali gases has drawn increasing attention. For 
example, in thermal equilibrium, both the spatial den- 
sity of the condensate atoms O and the distribution of 
number of particles in the condensate are modified due 
to the presence of the thermal atoms ||. Similarly, we 
must take into account the noncondensed atoms in or- 
der to explain time-dependent phenomena like the damp- 
ing and frequency shifts of collective modes |^,|| , or the 
evolution of the recently created vortices where the 
dissipation of the condensate motion is provided by the 
noncondensed atoms 0. 

The widely used Gross-Pitaevskii equation, the nonlin- 
ear Schrodinger equation for the condensate wave func- 
tion, is not able to describe these effects since it neglects 
the interaction between the condensate and the noncon- 
densed atoms H . One possibility to go beyond the Gross- 
Pitaevskii equation is to use Bogoliubov theory, which 
is a perturbative method valid for a small noncondensed 
fraction. The Bogoliubov method can be applied to ther- 
mal equilibrium but also to time-dependent situations: in 
a U(l) symmetry breaking point of view, to zeroth order 
one solves the time-dependent Gross-Pitaevskii equation 
for the condensate field vpo(r,t), to first order one lin- 
earizes the quantum field equations around the classical 
field ipo(r,t) to get the dynamics of noncondensed par- 
ticles, to second order one includes the back-action of 
noncondensed particles on the condensate. However, in 
this way one can predict only small corrections to the 
Gross-Pitaevskii equation. Additionally, if the number 
of noncondensed particles increases during the evolution 
of the system, the Bogoliubov approach is valid only for 
short times Another existing approach is the mean 
field Hartree-Fock-Bogoliubov approximation. This ap- 
proach is known however to present consistency problems 
and is still the object of research JlO| . 

In this paper we propose an alternative method to 



study the time evolution of Bose condensed gases per- 
turbed from an initial thermal equilibrium, based on the 
classical field approximation in the Wigner representa- 
tion, the so-called truncated Wigner approximation. The 
classical field approximation amounts to evolving a set 
of initially randomly distributed atomic fields with the 
usual Gross-Pitaevskii equation, with the crucial point 
that the field is now the whole matter field, not simply 
the condensate field. This approximation has been used 
already in the Glauber-P representation to study the for- 
mation of the condensate and it is expected to be correct 
for modes having a large occupation number Jll],[l2| ■ For 
such modes this approximation is known however to be 
most accurate in the Wigner representation as quantum 
noise is mimicked by classical noise in the initial state 
p3| . In addition in such representation, the approxi- 
mation is valid also for modes with a small occupation 
number as long as their evolution is well described by 
a quadratic approximation to the Hamiltonian, like in 
the Bogoliubov approach. The classical field approxima- 
tion is indeed exact in the Wigner representation for a 
quadratic Hamiltonian! The truncated Wigner approxi- 
mation is however expected to fail in describing the relax- 
ation of the gas to the correct thermal equilibrium, and 
one has to supplement the Gross-Pitaevskii equation by 
kinetic equations for the high energy modes Jll[ . 

The main result of this paper is the construction of 
an algorithm to approximately sample the Wigner dis- 
tribution at thermal equilibrium, using the number con- 
serving Bogoliubov theory This was the missing 
block for the implementation of the truncated Wigner ap- 
proach to study the dynamics of low temperature Bose 
condensed gases. We apply our method to the damping 
of the breathing mode of a ID condensate in a harmonic 
trap, which exemplifies the superiority of the truncated 
Wigner over the time-dependent Bogoliubov approach. 

We recall briefly the number conserving Bogoliubov 
approach that we use to describe the initial thermal equi- 
librium state of the gas. The total number of particles 
of the gas is supposed to be fixed and equal to N. The 
existence of the macroscopically populated mode <fi(r) 
motivates the splitting of the atomic field operator: 



(1) 



where (f> is the condensate wavefunction, d^ annihilates a 
particle in the mode <fi. The general idea is to eliminate 
the annihilation operator for the condensate. Unlike the 
usual approach, which replaces d^ by a fixed complex 



1 



number, we keep the operator nature of and write it 
in terms of the condensate phase operator and the 
number of condensate particles operator Nq as: 



a<f, = A<pN 



1/2 



(2) 



The modulus of can be re-expressed as a function of 
the number of noncondensed particles: 



N = N- / driP]_(r)iJj ± {r) 



(3) 



Since and A , are commuting unitary operators fli 



Id 



(4) 



commuting also with ip± , ip , we realize that the Hamil- 
tonian can be expressed as a function of the fields A(r ) = 
At^ip i_{r) and A^(r). In the Bogoliubov approximation 

we can then quadratize the Hamiltonian in terms of A: 
the terms linear in A vanish when <fi solves the time in- 
dependent Gross-Pitaevskii equation: 



H, 



- — A + C/(r) + 7V 5 |0| 2 - M 
2m 



= (5) 



where m is the atomic mass, U(f) is the trapping po- 
tential, g is the atomic coupling constant, and /j, is the 
chemical potential. The resulting quadratic form is the 
Bogoliubov Hamiltonian 



^ q uad(A) = J d?r^{k\-k)-C 



A 

At 



(6) 



where the differential operator C is given by: 



C 



QNg\<f>\ 2 Q 
-Q*Ng(r) 2 Q 



gp 



QNg<j) 2 Q* 
-H* gp - Q*Ng\ 



Here Q = Id — \4>)(<j>\ projects orthogonally to (j>. 

The initial A^-body density operator of the gas at tem- 
perature T is then approximated by 



a(t = 0) 



exp 



1 



k R T 



#quad(A) 



(8) 



The Wigner distribution associated to the A^-body den- 
sity operator a is a functional of a complex classical field 
ip(r). It is defined as the Fourier transform of a charac- 
teristic function X' 



X {l) = Tv\aeh^-^^ 



(9) 
(10) 



where J T> 2 j stands for the functional integral over the 
real and the imaginary part of the complex classical field 
7(r ) and the integrals in the exponent are over space. 



We now proceed with the analytical calculation of the 
Wigner distribution of the density operator (||). As we 
did with the atomic field operator, we split all the classi- 
cal fields into components parallel and orthogonal to the 
condensate mode </>, e.g. ip(f ) = ip\\<f>(f) + ip±('r)- 

It will prove convenient to use the commutation of 
with ip±_ to rewrite the characteristic function as 



x(l) 



!<{■ 



Til 



-h.c. J 7_ L A t J 4j ) -h.c. 



where {, } stands for the anticommutator and the brack- 
ets denote the expectation value in the density operator 
a . We then introduce an approximation for No in (|ll| ) 
replacing the operator N in Eq. (||) by its actual value N 
in the gas: 



N ~ N - 



A+A. 



(12) 



As a consequence the phase operators A^,A^ commute 
with all the operators appearing in X- We now calculate 
the trace in (|ll]) using a simultaneous eigenbasis of the 
phase operator A,/,, with eigenvalue e l6 , and of the Bo- 
goliubov Hamiltonian whose eigenvalues depend on the 
occupation numbers of elementary excitations but not on 
the phase 9. The form obtained for W clearly preserves 
the U{\) symmetry: 



W(tp) 



2ji 



(13) 



Wo is the 9 — component of the Wigner distribution. 
By integrating over the real and the imaginary parts of 
7|| and by using the cyclic properties of the trace we get 



(7) WoW) = S^j) Wig A \\ [styf - # V V} 



W ± ) (14) 



Ft I 

where ip« ' are the real and imaginary parts of tpu • The 



expression Wig^ 



/(A) 



II 

denotes the Wigner transform 



with respect to A of any function / of A. It is obtained 
from (||) and ( |l0| ) by the substitutions: 



/(A) $ -> A ip-*tp± 7-*7±. 



(15) 



It is convenient to rewrite Wo(ip) as a quasiprobability 
distribution for Nq = ipn?P\\ and ip±: 



P(NoM = Wig A 



^{6(N Q -N Q ),&} 



(i> ± ). (16) 



A first insig ht into @ can be obtained by looking at its 
marginals. By integrating over Nq: 



p(V>j 



dN P(N M = Wig A (<x)(V>j 



(17) 



2 



we obtain the Wigner distribution for noncondensed 
modes, which is a Gaussian IflSfl: 



P(tp±) oc exp ^ - / <fr (xj;* ± , tp x ) ■ M 



M 



\ ? | tanh — — — 
-1 J 2k B T 



(18) 



(19) 



where C is the Bogoliubov operator (|t]). The other 
marginal is obtained by integrating ( |l6| ) over ip±: 

P(N ) = J V 2 t/j ± P(N , 1> x ) = Tr[S(N - N ) &]. (20) 

P(Nq) is then simply the probability distribution of the 
number of particles in the condensate Jl{| . 

Let us now discuss the limits of validity of ([l6]). The 
spectrum of Nq contains integer numbers only so that 
S(Nq — No) in Eq.([l6|) vanishes when Nq is not integer. 
As a function of No, P(Na,ip±) is therefore a "comb" 
of delta functions centered on integer values. The exact 
Wigner distribution should be instead a smooth function 
of its variables. The comb-artifact comes from the ap- 
proximation @ that we have made on the condensate 
mode. Our description is sensible if there are many peaks 
within the width of the distribution of No for a given ip± . 
The variance of Ao for a given ip± , that is the conditional 
variance of No, should then satisfy: 



Var(A |V>_L) > 1. 



(21) 



We estimate this variance for a "typical" ip± by taking 
the mean of Var(iVo|^j_) over the probability distribution 
of ip± given in (|l8|): 



(Var(iVolV'x))^ = ~ Tr± [id - M 



(22) 



where the trace is over the Bogoliubov modes. For the 
modes with energy much larger than ksT one has M 2 ~ 
Id, while for the modes with energy lower or equal to 
ksT one has M 2 <C Id, so that: 



(Var(iVo|Vx))vu 



1 



K 



th 



(23) 



where A/th is the number of thermally populated modes 
which we require to be large. In practice ksT should be 
high enough: for a condensate in an isotropic harmonic 
trap of frequency uj, we should have fc^T 3> hu). 



The next point is to sample the Wigner function (16) 



We proceed in two steps: (i) First we sample tfjj_ irre- 
spectively to the value N according to the Gaussian dis- 
tribution (|l8|) with a stochastic field method Jl^,|2(J . (ii) 
We have to sample the conditional distribution of No for 
a fixed ip±: 



It turns out that the width in No of the distribution 
P(No\4>±) is much narrower than the width of P(Nq) 
when the condition ( pl| ) is satisfied. In fact as soon as 
a sufficient number of modes is populated, the following 
inequality holds involving the number of modes A/th with 
energy lower or equal to ksT, the average of the num- 
ber of particles out of the condensate (N ~ No), and its 
variance Var(A — No): 

Mh <(N- No) < Vax(N - N Q ) = Var(A ) . (25) 

We then replace ( pij ) with a delta distribution centered 
on the mean value of Ao for a fixed ip±: 

Mean(Aol^) = C - \ J [W - M 2 } 

where C is a constant plj |. 

Finally, for each randomly generated ip± and Ao we 
are now able to form the total field ip as: 



^ = A o 1/2 (0+^] -:■ r 



(26) 



In @ the function 4>^ /N is a correction to the con- 
densate wavefunction induced by the back-action of non- 
condensed particles onto the condensate. This correction 
has to be included since it gives a contribution to observ- 
ables on the same order as ip±_ pq |. We calculate (f)^ /N 
using the procedure of Q . 

In the classical field approximation each field if) evolves 
according to the time-dependent Gross-Pitaevskii equa- 
tion Onl: 



ih d t ifj 



2m 



A + U(r,t)+g\iP\ 2 



(27) 



P(N \tp±) = P(N ,i> ± )/P(ip ± ). 



(24) 



The average of an observable at time t will be obtained 
by averaging over many stochastic classical fields i/j(t). 

We now give an illustration of our method. In Fig- 
ure |] we show density oscillations in the center of a ID 
Bose condensed gas induced by an abrupt change of the 
harmonic trap frequency uj — > 0.8a;. While the Gross- 
Pitaevskii equation for the condensate field predicts un- 
damped oscillations, we get damping at finite tempera- 
ture when we average over stochastic realizations in the 
Wigner method J2^]. The time-dependent Bogoliubov 
theory in [ fl6| correctly predicts the damping at short 
times but gives unphysically large oscillations at longer 
times. We understand this failure as due to the growth in 
time of the number of noncondensed particles and of the 
condensate wavefunction correction <fi( 2 '. This growth 
invalidates the linearized treatment at the basis of the 
time-dependent Bogoliubov approach fl , at a relatively 
short time here because of the small number of particles. 
On the contrary, in the truncated Wigner approach, no 
linearized treatment is performed and the whole field %j) 
is evolved according to the fully non-linear equation (p7|). 
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FIG. 1. Damped oscillations of the mean density in the 
center of a ID condensed cloud in a harmonic trap after an 
abrupt change of the trap frequency u) — > 0.8a;. fcsT = 307kj, 
fj, — S.lhui and N = 10 3 . Solid line: Wigner simulation. 
Dotted curve: Bogoliubov theory, p is in units of (mui/H) 1 ' 2 . 

In Figure 1 a spatial grid of 64 point has been used 
for discretization in the numerical implementation of our 
algorithm. We have tested that the results are unchanged 
by doubling the number of points in the grid [ p3| . 

In conclusion we have developed an efficient algorithm 
for approximate stochastic sampling of the Wigner rep- 
resentation of the TV-body density matrix of a Bose con- 
densed gas in thermal equilibrium. The initial distribu- 
tion of atomic fields ip given by this sampling can then be 
evolved in the classical field (truncated Wigner) approxi- 
mation to study the response of the gas to an excitation. 
As far as time-independent situations are concerned, our 
approach is equivalent to the number conserving Bogoli- 
ubov approach. In the time dependent case however un- 
like the Bogoliubov approach, the classical field evolution 
is fully nonlinear both for the condensate and the non- 
condensed particles, which extends the applicability of 
the method to longer evolution times. 
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